rm(list = ls())

library(bcp)
library(strucchange)
library(foreign)

bg <- read.dta("updateddividedsenate.dta")

pctdiv.ts <- ts(bg$percentdivided)

pctdiv <- bg$percentdivided

## compute breakpoints--note that we are regressing time series of
## coalition size on a constant--we could include covariates 

n <- length(pctdiv.ts)

## priors:

## f(p) 1 /p_0, O < p < p_0,

## f(w) =I/w_0, O < w < wO,

bcp.pctdiv <- bcp(pctdiv, p0 = .2, w0 = .2, mcmc=50000)

yrlabs <- seq(1788,2004,2)

yrlabs4 <- seq(1788,2004,4)

bcp.results <- cbind(yrlabs,bcp.pctdiv$posterior.prob,bcp.pctdiv$posterior.mean)

colnames(bcp.results) <- c("year","Post. Prob.","X1")

bcp.results

plot.bcp(bcp.pctdiv)

pdf (paste(file="pctdiv_bcp_bp_plot_nolags_h3.pdf"))

bp.pctdiv <- breakpoints(pctdiv ~ 1, h = 3, breaks=15)

## note that h (minimal segment length), must be > than number of regressors

ci.pctdiv <- confint(bp.pctdiv, level = 0.90)

summary(bp.pctdiv)

ci.pctdiv2 <- confint(bp.pctdiv, level = 0.90, het.reg=T, het.err=F)

x <- bcp.pctdiv

posterior.prob <- x$posterior.prob

posterior.prob[length(posterior.prob)] <- 0
		
op<-par(mfrow=c(2,1),col.lab="black",col.main="black")

op2 <- par(mar=c(0,4,4,2),xaxt="n", cex.axis=0.75)

plot(1:length(x$data), x$data, col="grey", pch=20, xlab="", ylab="Posterior Mean", main="Posterior Means and Probabilities of a Change")

t <- seq(1,length(pctdiv),1)

t4 <- seq(1,length(pctdiv),2)

lines(x$posterior.mean, lwd=2)

ci.pctdiv$confint

ci.pctdiv2$confint

for (i in 1:nrow(ci.pctdiv$confint)) {

  abline(v=ci.pctdiv2$confint[i,1],col=2,lty=3)
  abline(v=ci.pctdiv2$confint[i,3],col=2,lty=3)
  abline(v=ci.pctdiv2$confint[i,2],col=1,lty=2)
  
}

op3 <- par(mar=c(5,4,0,2), xaxt="s", cex.axis=0.75)

plot(1:length(x$posterior.mean), posterior.prob, type="l",
     ylim=c(0,max(posterior.prob)), xlab="Year",
     ylab="Posterior Probability", main="",xaxt="n")

par(cex.axis=.7)

axis(1, at=t4,labels=yrlabs4,mgp=c(2,.7,0),las=2)

for (i in 1:nrow(ci.pctdiv2$confint)) {

  abline(v=ci.pctdiv2$confint[i,1],col=2,lty=3)
  abline(v=ci.pctdiv2$confint[i,3],col=2,lty=3)
  abline(v=ci.pctdiv2$confint[i,2],col=1,lty=2)
  
}

cbind(seq(1,length(bg$year),1),bg$year)
